Biometrics 1 (Bio360) Final Test

A case study based on

Ural, B.B., Caron, D.P., Dogra, P. et al. Inhaled particulate accumulation with age impairs immune function and architecture in human lung lymph nodes. Nat Med 28, 2622–2632 (2022). https://doi.org/10.1038/s41591-022-02073-x

library(readr)
library(ggplot2)
library(purrr)
library(dplyr)

Study Part 1

Lung lymphs nodes (LLNs) accumulate carbon-pollutant particulates and turn black with age, unlike mesenteric lymph nodes (MLNs).

We make two linear regression models, one for LNNs and one for MLNs.

x = Age y = Percent Lymph Node showing particulate accumulation

A.da <- read_csv("AgeLNdata.csv")

Rows: 35 Columns: 3── Column specification ─────────────────────────────────────────────────────
Delimiter: ","
chr (1): Type
dbl (2): Age, Percent
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

ggplot(A.da, aes(x=Age, y=Percent, shape=Type, color=Type)) +
  geom_point() +
  geom_smooth(method=lm, se=FALSE, fullrange=TRUE) +
  labs(title= "Age and Lympth Node Particulate Accumulation", y= "% Particulate Accumulation") +
  theme_bw() + theme(text = element_text(size =13))  

ggsave("AgeVsLNs.pdf")

Saving 6.12 x 3.79 in image

ggsave("AgeVsLNs.png")

Saving 6.12 x 3.79 in image

Subset data

LLN<- A.da[A.da$Type=="LLN",]
MLN <-A.da[A.da$Type=="MLN",]

mod.LLN <-lm(Age~Percent, data=LLN)
summary(mod.LLN)

Call:
lm(formula = Age ~ Percent, data = LLN)

Residuals:
    Min      1Q  Median      3Q     Max 
-23.885 -16.533  -3.897  10.950  37.501 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  34.8598     5.8503   5.959 1.56e-05 ***
Percent       1.2404     0.4047   3.065    0.007 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 19.74 on 17 degrees of freedom
Multiple R-squared:  0.356, Adjusted R-squared:  0.3181 
F-statistic: 9.397 on 1 and 17 DF,  p-value: 0.007004

Check residuals

par(mfrow=c(2,2))
plot(mod.LLN)
par(mfrow=c(1,1))

mod.MLN <-lm(Age~Percent, data=MLN)
summary(mod.MLN)

Call:
lm(formula = Age ~ Percent, data = MLN)

Residuals:
    Min      1Q  Median      3Q     Max 
-32.523 -17.273  -3.015  19.038  37.165 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   43.523      7.802   5.579  6.8e-05 ***
Percent       22.133     44.348   0.499    0.625    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 24.87 on 14 degrees of freedom
Multiple R-squared:  0.01748,   Adjusted R-squared:  -0.0527 
F-statistic: 0.2491 on 1 and 14 DF,  p-value: 0.6255

par(mfrow=c(2,2))
plot(mod.MLN)
par(mfrow=c(1,1))

Part 2

The CD68+CD169− subset of LN macrophages contain the particulates.

Use an ANOVA to test if one type of macrophage disproportionately accumulates these pollutant particulates.

Cell Count is # of cells that contain particulates/mm2

da2 <- read_csv("part2data.csv")

Rows: 54 Columns: 3── Column specification ─────────────────────────────────────────────────────
Delimiter: ","
chr (2): Age, Macrophage
dbl (1): CellCount
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

da2$Macrophage <- as.factor(da2$Macrophage)

ggplot(da2, aes(x=Macrophage, y=CellCount)) + 
    geom_boxplot() +
      theme_bw() + theme(text = element_text(size = 14))  

ggsave("ANOVAboxplot.pdf")

Saving 6.12 x 3.79 in image

ggsave("ANOVAboxplot.png")

Saving 6.12 x 3.79 in image

Using

M.stats <- da2 %>% 
  split(.$Macrophage) %>% 
  map(summary)
M.stats

$MP1.
     Age            Macrophage   CellCount      
 Length:18          MP1.:18    Min.   :  37.87  
 Class :character   MP2.: 0    1st Qu.: 329.57  
 Mode  :character   MP3.: 0    Median : 511.73  
                               Mean   : 624.35  
                               3rd Qu.: 793.79  
                               Max.   :1945.99  

$MP2.
     Age            Macrophage   CellCount     
 Length:18          MP1.: 0    Min.   : 11.68  
 Class :character   MP2.:18    1st Qu.: 45.54  
 Mode  :character   MP3.: 0    Median : 62.65  
                               Mean   :107.80  
                               3rd Qu.:158.28  
                               Max.   :391.38  

$MP3.
     Age            Macrophage   CellCount     
 Length:18          MP1.: 0    Min.   :  0.32  
 Class :character   MP2.: 0    1st Qu.: 10.32  
 Mode  :character   MP3.:18    Median : 35.91  
                               Mean   : 45.05  
                               3rd Qu.: 79.39  
                               Max.   :117.55  

tapply(da2$CellCount, da2$Macrophage, c(summary, mean))    # Summary by group using tapply

Error in match.fun(FUN) : 
  'c(summary, mean)' is not a function, character or symbol

da2 %>% 
  select(-Age) %>% 
  group_by(Macrophage) %>% 
  summarise(mean=mean(CellCount), sd=sd(CellCount), n=n())

NA
NA

AnovaMod <- aov(CellCount ~ Macrophage, data=da2)
summary(AnovaMod)

            Df  Sum Sq Mean Sq F value   Pr(>F)    
Macrophage   2 3638148 1819074   24.64 3.26e-08 ***
Residuals   51 3765551   73834                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

lsr::etaSquared(AnovaMod)

             eta.sq eta.sq.part
Macrophage 0.491396    0.491396

tukey.res<- TukeyHSD(AnovaMod)
tukey.res

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = CellCount ~ Macrophage, data = da2)

$Macrophage
                diff       lwr       upr     p adj
MP2.-MP1. -516.54944 -735.1954 -297.9035 0.0000018
MP3.-MP1. -579.30667 -797.9526 -360.6607 0.0000001
MP3.-MP2.  -62.75722 -281.4032  155.8888 0.7687074

#plotTukeyHSD(tukey.res)

plot(tukey.res)

par(mfrow=c(2,2))
plot(AnovaMod)
par(mfrow=c(1,1))

Regression Models

data <- read_csv("regressiondata.csv")

Rows: 54 Columns: 3── Column specification ─────────────────────────────────────────────────────
Delimiter: ","
chr (1): Macrophage
dbl (2): Age, Arginine_IL6
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

ggplot(data, aes(x=Age, y=Arginine_IL6, color=Macrophage)) +
  geom_point(size=2.5) +
  geom_smooth(method=lm, se=FALSE, fullrange=TRUE) +
  labs(title= "Regression Model Matching") +
  theme_bw() +theme(text = element_text(size = 12))  


ggsave("RegModelMatch.pdf")

Saving 6.12 x 3.79 in image

ggsave("RegModelMatch.png")

Saving 6.12 x 3.79 in image

a <- data[data$Macrophage == "CD3. w/particulates",]

mod.a <- lm(Arginine_IL6 ~ Age, data=a)
summary(mod.a)

Call:
lm(formula = Arginine_IL6 ~ Age, data = a)

Residuals:
    Min      1Q  Median      3Q     Max 
-16.652 -11.293  -4.426   0.565  39.175 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  -7.0466    20.2848  -0.347    0.735
Age           0.3374     0.3212   1.051    0.316

Residual standard error: 18.28 on 11 degrees of freedom
Multiple R-squared:  0.09119,   Adjusted R-squared:  0.008575 
F-statistic: 1.104 on 1 and 11 DF,  p-value: 0.316

b <- data[data$Macrophage == "CD3. w/o particulates",]

mod.b <- lm(Arginine_IL6 ~ Age, data=b)
summary(mod.b)

Call:
lm(formula = Arginine_IL6 ~ Age, data = b)

Residuals:
    Min      1Q  Median      3Q     Max 
-17.064  -8.169  -2.702   7.953  24.112 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) 18.31882   14.68600   1.247    0.238
Age          0.01869    0.23253   0.080    0.937

Residual standard error: 13.23 on 11 degrees of freedom
Multiple R-squared:  0.0005867, Adjusted R-squared:  -0.09027 
F-statistic: 0.006458 on 1 and 11 DF,  p-value: 0.9374

c <- data[data$Macrophage == "CD2. w/particulates",]

mod.c <- lm(Arginine_IL6 ~ Age, data=c)
summary(mod.c)

Call:
lm(formula = Arginine_IL6 ~ Age, data = c)

Residuals:
    Min      1Q  Median      3Q     Max 
-16.894 -11.451  -6.949   4.206  42.530 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  24.7018    14.4540   1.709    0.113
Age          -0.1421     0.2638  -0.539    0.600

Residual standard error: 17.27 on 12 degrees of freedom
Multiple R-squared:  0.02362,   Adjusted R-squared:  -0.05774 
F-statistic: 0.2903 on 1 and 12 DF,  p-value: 0.5999

d <- data[data$Macrophage == "CD2. w/o particulates",]

mod.d <- lm(Arginine_IL6 ~ Age, data=d)
summary(mod.d)

Call:
lm(formula = Arginine_IL6 ~ Age, data = d)

Residuals:
    Min      1Q  Median      3Q     Max 
-34.117 -20.263   1.662  16.014  40.096 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  44.9565    20.0834   2.238   0.0449 *
Age          -0.2004     0.3665  -0.547   0.5945  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 24 on 12 degrees of freedom
Multiple R-squared:  0.02432,   Adjusted R-squared:  -0.05699 
F-statistic: 0.2991 on 1 and 12 DF,  p-value: 0.5945

Part 3

Particulate Accumulation impairs immune function

“Expression of the key activation markers CD80 and CD86 and the phagocytic marker CD36 decreased with age specifically in CD68+CD169− macrophages in LLNs but not in MLNs. In contrast, CD209 expression was not altered significantly with age in any macrophage subset at either site (Fig. 3e and Extended Data Fig. 4b). These results show that while the frequencies of LN macrophage subsets are largely maintained as age progresses, the expression of functional markers specifically in the CD68+CD169− subset within LLNs decreases with age, suggesting that particulates may have specific effects on macrophage function.”

part3data <- read_csv("part3data.csv")

Rows: 120 Columns: 5── Column specification ─────────────────────────────────────────────────────
Delimiter: ","
chr (3): Marker, MacrophageType, LNType
dbl (2): Age, Percent
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

part3data$Label <- "X"
part3data[42,6] <- "P3"
part3data[52,6] <- "P4"
part3data[41,6] <- "P1"
part3data[33,6] <- "P2"
part3data$Label <- as.factor(part3data$Label)

part3dataLLN <- part3data[part3data$LNType == "LLN",]
part3dataMLN <- part3data[part3data$LNType == "MLN",]

geom_text(data=subset(mtcars, wt > 4 | mpg > 25), aes(wt,mpg,label=name))

ggplot(part3dataLLN, aes(x=Age, y=Percent, shape=Marker, color=Marker, label= Label)) +
  geom_point() +
  geom_smooth(method=lm, se=FALSE, fullrange=TRUE) +
  labs(title= "LLN immunity markers and age", y= "% Marker") +
  theme_bw()+
  geom_text(data = subset(part3dataLLN,
                          Label == "P1"| Label == "P2"| Label == "P3" | Label=="P4"), aes(label=Label), hjust=0, vjust=1)

ggsave("LNNmarkers.pdf")

Saving 7.29 x 4.51 in image

ggsave("LNNmarkers.png")

Saving 7.29 x 4.51 in image

da36 <- part3dataLLN[part3dataLLN$Marker == "CD36",]

modCD36LLN <- lm(Percent ~ Age, data=da36)
summary(modCD36LLN)

Call:
lm(formula = Percent ~ Age, data = da36)

Residuals:
    Min      1Q  Median      3Q     Max 
-18.098 -10.660  -1.013   3.273  29.210 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  62.0207     7.6238   8.135 1.93e-07 ***
Age          -0.6278     0.1342  -4.678 0.000187 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 14.13 on 18 degrees of freedom
Multiple R-squared:  0.5487,    Adjusted R-squared:  0.5236 
F-statistic: 21.88 on 1 and 18 DF,  p-value: 0.0001873

Check residuals

par(mfrow=c(2,2))
plot(modCD36LLN)
par(mfrow=c(1,1))

da209 <- part3dataLLN[part3dataLLN$Marker == "CD209",]

mod209LLN<- lm(Percent ~ Age, data=da209)
summary(mod209LLN)

Call:
lm(formula = Percent ~ Age, data = da209)

Residuals:
    Min      1Q  Median      3Q     Max 
-20.357 -10.186  -5.598   8.507  39.633 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   4.8753     9.1474   0.533     0.60
Age           0.2466     0.1605   1.537     0.14

Residual standard error: 16.96 on 20 degrees of freedom
Multiple R-squared:  0.1057,    Adjusted R-squared:  0.06094 
F-statistic: 2.363 on 1 and 20 DF,  p-value: 0.1399

par(mfrow=c(2,2))
plot(mod209LLN)
par(mfrow=c(1,1))

da80 <- part3dataLLN[part3dataLLN$Marker == "CD80/86",]

modCD80LLN <- lm(Percent ~ Age, data=da80)
summary(modCD80LLN)

Call:
lm(formula = Percent ~ Age, data = da80)

Residuals:
    Min      1Q  Median      3Q     Max 
-15.379  -6.078  -0.011   2.058  42.581 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  26.9513     6.7961   3.966 0.000829 ***
Age          -0.3012     0.1194  -2.522 0.020746 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 12.6 on 19 degrees of freedom
Multiple R-squared:  0.2508,    Adjusted R-squared:  0.2114 
F-statistic: 6.361 on 1 and 19 DF,  p-value: 0.02075

par(mfrow=c(2,2))
plot(modCD80LLN)
par(mfrow=c(1,1))

da209b <- da209[-c(1, 11),]

mod209LLNb<- lm(Percent ~ Age, data=da209b)
summary(mod209LLNb)

Call:
lm(formula = Percent ~ Age, data = da209b)

Residuals:
    Min      1Q  Median      3Q     Max 
-12.403  -6.970  -4.124   9.395  24.411 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept) 11.59479    5.80831   1.996   0.0613 .
Age          0.03386    0.10644   0.318   0.7541  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 10.35 on 18 degrees of freedom
Multiple R-squared:  0.00559,   Adjusted R-squared:  -0.04965 
F-statistic: 0.1012 on 1 and 18 DF,  p-value: 0.7541

ggplot(da209b, aes(x=Age, y=Percent, shape=Marker, color=Marker)) +
  geom_point() +
  geom_smooth(method=lm, se=FALSE, fullrange=TRUE) +
  labs(title= "LLN immunity markers and age", y= "% Marker") +
  theme_bw()

modcomLLN<- lm(Percent ~ Age, data=da36.8086)
summary(modcomLLN)

Call:
lm(formula = Percent ~ Age, data = da36.8086)

Residuals:
    Min      1Q  Median      3Q     Max 
-29.701 -10.636  -3.895   6.119  37.962 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  44.4247     6.2797   7.074 1.69e-08 ***
Age          -0.4670     0.1104  -4.229 0.000137 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 16.46 on 39 degrees of freedom
Multiple R-squared:  0.3144,    Adjusted R-squared:  0.2968 
F-statistic: 17.88 on 1 and 39 DF,  p-value: 0.0001373

ggplot(part3dataMLN, aes(x=Age, y=Percent, shape=Marker, color=Marker)) +
  geom_point() +
  geom_smooth(method=lm, se=FALSE, fullrange=TRUE) +
  labs(title= "MLN immunity markers and age", y= "% Marker") +
  theme_bw()

ggsave("MLNmarkers.pdf")

Saving 7.29 x 4.51 in image

ggsave("MLNmarkers.png")

Saving 7.29 x 4.51 in image

da36M <- part3dataMLN[part3dataMLN$Marker == "CD36",]

modCD36MLN <- lm(Percent ~ Age, data=da36M)
summary(modCD36MLN)

Call:
lm(formula = Percent ~ Age, data = da36M)

Residuals:
    Min      1Q  Median      3Q     Max 
-25.093  -5.591   0.518   2.512  31.508 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  47.6271     9.2026   5.175  7.6e-05 ***
Age          -0.0456     0.2000  -0.228    0.822    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 13.01 on 17 degrees of freedom
Multiple R-squared:  0.003049,  Adjusted R-squared:  -0.05559 
F-statistic: 0.052 on 1 and 17 DF,  p-value: 0.8223

par(mfrow=c(2,2))
plot(modCD36MLN)
par(mfrow=c(1,1))

da209M <- part3dataMLN[part3dataMLN$Marker == "CD209",]

mod209MLN <- lm(Percent ~ Age, data=da209M)
summary(mod209MLN)

Call:
lm(formula = Percent ~ Age, data = da209M)

Residuals:
    Min      1Q  Median      3Q     Max 
-19.564  -7.830  -2.559   7.182  29.784 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  30.7766     8.3685   3.678  0.00187 **
Age          -0.1619     0.1819  -0.890  0.38582   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 11.83 on 17 degrees of freedom
Multiple R-squared:  0.04453,   Adjusted R-squared:  -0.01167 
F-statistic: 0.7923 on 1 and 17 DF,  p-value: 0.3858

par(mfrow=c(2,2))
plot(mod209MLN)
par(mfrow=c(1,1))

da80M <- part3dataMLN[part3dataMLN$Marker == "CD80/86",]

mod80MLN <- lm(Percent ~ Age, data=da80M)
summary(mod80MLN)

Call:
lm(formula = Percent ~ Age, data = da80M)

Residuals:
    Min      1Q  Median      3Q     Max 
-16.358  -8.367  -3.173   2.842  36.045 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept) 36.28337    9.58759   3.784  0.00148 **
Age          0.09941    0.20836   0.477  0.63934   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 13.56 on 17 degrees of freedom
Multiple R-squared:  0.01321,   Adjusted R-squared:  -0.04483 
F-statistic: 0.2277 on 1 and 17 DF,  p-value: 0.6393

par(mfrow=c(2,2))
plot(mod80MLN)
par(mfrow=c(1,1))